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Kinetic equations, which explicitly take into account the 
branching nature of sandpile avalanches, are derived. The 
dynamics of the sandpile model is described by the generat- 
ing functions of a branching process. Having used the results 
obtained the renormalization group approach to the critical 
behavior of the sandpile model is generalized in order to cal- 
culate both critical exponents and height probabilities. 

PACS numbers: 05.40. +j, 05.60.+W, 02.50.-r, 64.60. Ak 



As it was realized a few years ago, almost all fractal 
structures we can see in nature are "self-organized" in a 
sense that they spontaneously grow, obeying rather sim- 
ple dynamical rules 0. The concept of Self-Organized 
Criticality (SOC) has recently been introduced by Bak 
et al. |H to emphasize the fact that the production of 
self-similar structures is not the only possible result of 
the irreversible dynamics of an extended system. They 
have shown that such a system can also evolve stochasti- 
cally into a certain critical state similar to that of second 
order phase transition. It lacks therein any characteristic 
length as well as timescale and obeys power-law distri- 
butions. The critical state is independent of the initial 
configuration of the system and, unlike ordinary critical 
phenomena, no fine tuning of control parameters is nec- 
essary to arrive at this state. 

To illustrate the basic ideas of SOC, Bak et al. used a 
cellular automaton now commonly known as "sandpile" 
because of the crude analogy between its dynamical rules 
and the way sand topples when building a real sand pile. 
The formulation of this model is given in terms of integer 
height variables Zi at each site of a square lattice C. Par- 
ticles are added randomly and the addition of a particle 
increases the height at that site by one. If this height ex- 
ceeds the critical value z c = 4, then the site topples, and 
on toppling its height decreases by 4 and the heights at 
each of its nearest neighbors increases by 1. They may 
become unstable in their turn and the dynamical pro- 
cess continues. Open boundary conditions are usually 
assumed, so that particles can leave the system. 

Dhar H has shown that this model is exactly solvable 
due to an Abelian group structure hidden in its dynam- 
ics. The fact that almost all characteristics of the Abelian 
sandpile model are calculable analytically is the reason 
for considering this model as a perfect proving ground for 
various approaches to SOC. Among others, the Renor- 
malization Group (RG) approach proposed recently by 
Pietronero et al. |^,^) seems to be the most promising one 



because it explicitly describes the self-similar dynamics 
of the SOC models at the critical state. 

The purpose of this Letter is to generalize this renor- 
malization scheme by exploiting the analogy between 
large scale dynamics of the sandpile model and chain 
reactions of a special kind. At first, we should define 
coarse grained variables proper for describing the behav- 
ior of the system on a sublattice Cb which consists of the 
cells of size b on the initial lattice C. To keep the connec- 
tion with the original formulation of the sandpile model, 
we will characterize the static properties of a cell by four 
quantities 

n (b) = (n A , n B , n c , «d), n A + n B + n c + = 1, (1) 

which are nothing but the probabilities for a cell to be- 
have like a site on the initial lattice with a height 1, 2, 3 
or 4 respectively in the coarse grained dynamics, i.e. the 
addition of a "coarse grained particle" to the cell trans- 
forms it to the next one in the alphabet. For example, 
the cell B characterized by the vector (0,1,0,0) will be 
transformed to the cell C with the vector (0,0,1,0). The 
last variable ?id is the probability for the cell to behave 
like a critical one in a sense that the addition of a " coarse 
grained particle" to the cell induces relaxations into some 
neighboring cells or, in other words, subrelaxation pro- 
cesses on a minimal scale span the cell and transfer energy 
to some of its neighbors. 

According to Pietronero et al. 0] , independently of the 
dynamics of the model at the minimal scale such a relax- 
ation process leads to four possible situations for a coarse 
grained cell of size b. Namely, after the relaxation of the 
critical cell "coarse grained particles" can be transferred 
to one, two, three or four neighboring cells with the prob- 
abilities 

P (fc) = (Pl,P2,P3,Pi), Pi +P2 +P3 +P4 = 1. (2) 

The probabilities pi are normalized sums of all the pro- 
cesses that ignite the corresponding number of neighbor- 
ing cells, independently of their position. The distribu- 
tion of particles after toppling at the minimal scale is 
characterized by the vector pW = (0, 0, 0, 1). 

While the first set of the coarse grained variables 
describes the height configuration at the scale b, the sec- 
ond set p*- 6 ) defines the toppling rules of the model and, 
in some sense, characterizes the phase space for the re- 
laxation dynamics at the same scale. 

In this framework the coarse grained dynamics of 
the sandpile model can be represented as the following 
branching process on the sublattice Cb 
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Here <p and </5 denote the "coarse grained particles" ob- 
tained by the cell and the particles transferred to the 
neighboring cells, respectively. 

These processes can formally be reinterpreted as an 
irreversible chemical reaction which takes place at each 
cell of the sublattice Now the coarse grained vari- 
ables riA , riB , «c , ™d and n v denote the concentrations 
of the respective species A, B, C, D, and ip. Following 
to standard prescriptions of the chemical physics we can 
write kinetic equations corresponding to this scheme of 
chemical reactions 
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where p — p\ + 2p2 + 3p% + Ap± is equal to the average 
number of particles leaving the cell on toppling and r 
is the position vector of the cell in the 2D space. The 
noise term r]{Y,t), being non-negative, mimics the ran- 
dom addition of particles to the system. The diffusion 
term V 2 (n lp nu) describes the transfer of particles into 
the neighboring cells, and the diffusion coefficient v for 
the discrete Laplacian on the square lattice is equal to 
1/4. 

The only mobile specie in this scheme of reactions is 
ip and it is the field n v which describes the dynamics of 
avalanches. When it is equal to zero, all toppling pro- 
cesses die. Then, due to the noise term n(r,t), particles 
are added randomly into the system initiating a branch- 
ing process directed to the open boundary of the system. 
This process mutates species in the cells it has visited 
and topples the critical ones. Finally, the system will 
reach the steady state where the probability that the ac- 
tivity will die is on average balanced by the probability 
that the activity will branch. Thus, the chain reaction 
maintains this stationary state and all further avalanches 
cannot change the concentrations of species A, B, C, and 
D. Therefore, the steady state is characterized by the 
conditions that 



n-A = h B = h c = hu = 



(5) 



and Eqs. (4a-4d|) lead to the following relationships be- 
tween concentrations of species n( fc ) at the stationary 
state and branching probabilities pW 



n *A =Pi/P; 

n B = (P3+P4)/P, 

n* c = {P2+P3+Pi)/P, 

n h = (Pi +P2+P3+ PA)Ip = l/P 
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The relation (^) between the probability and branch- 
ing probabilities has already appeared in the paper 
Q. It was derived there from the assumption that at 
the stationary state the flow of particles in a cell was on 
average balanced by the flow of particles out of the cell. 

If we neglect the fluctuations of n( b ) at the steady state, 
Eq. (|4e|), describing the propagation of an active process, 
becomes simply the diffusion equation 



(7) 



and coincides with that for the flow of particles at the 
critical state obtained by Zhang ||. 

To describe in detail the branching process underly- 
ing the large scale behavior of the sandpile model, let us 
consider the following generating function: 
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where symbols N, E, S and W correspond to the north, 
east, south and west directions on the square lattice, re- 
spectively. The coefficient of each term of this polyno- 
mial, say P2NE/6, is equal to the probability that after 
relaxation of a critical cell particles will go only north- 
ward and eastward. 

It is easy to check directly that this function has the 
following properties: (a) if any letter, for example N, 
is replaced by zero, the function tr(0, E, S, w) counts all 
the processes that do not send a particle northward; (b) 
on the contrary, the difference cr(l, E, S, w) — cr(0, E, S, w) 
corresponds to the sum over all the processes that defi- 
nitely send a particle to the north and, possibly, to some 
other directions; (c) this generating function is normal- 
ized so that (7(1,1,1,1) = 1 due to the normalization 
condition on branching probabilities, Eq. (||); (d) if the 
particle that has been sent northward after the relax- 
ation of a critical cell ignites the northern neighboring 
cell, the corresponding letter N should be replaced by an- 
other cr-function and each term of the generating function 
cr(cr(N, E, S, w), E, S, w) - cr(0, E, S, w) corresponds to the 
process that consists of the two successive toppiings. 

Using these simple rules we can write the generating 
function for any chain of relaxation processes on the sub- 
lattice Lb- 

Now, following to the general ideas of the paper [Q, 
we are ready to define a renormahzation transformation 
for the relaxation dynamics of the sandpile model. The 
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standard real space RG approach consists of considering 
a block 2 x 2 of cells of the lattice Lb to be a single cell 
of the size 2b on the lattice £26- Thus, given branch- 
ing probabilities we want to find an analogous set 
of probabilities p( 2fc ) on the lattice 11,2b- To this end, we 
should count up all the possible toppling processes that 
span the starting block of cells and transfer particles to 
some neighboring blocks. The spanning rule implies that 
we have to consider only those connected chains of pro- 
cesses that span the block from left to right or from top to 
bottom neglecting the processes not extending over the 
resulting scale 26. On Fig. [l] all the different types of top- 
pling processes that satisfy this spanning condition are 
shown. The blocks that have no, or have only one criti- 
cal cell inside are not included in the renormalization of 
the dynamics because they do not lead to a relaxation 
process that spans the cell of size 2b. 

We amplify the set of toppling processes considered 
in the RG scheme of Pietronero et al. [Q by taking into 
account also the processes, Fig. [j] (c), of relaxation of 
the cell C provided it has got two particles during the 
relaxation of the block. These additional processes seem 
to be very important for the relaxation dynamics of the 
sandpile model, and taking account of them improves the 
accuracy of the RG approach considerably. 
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FIG. 1. We show the four different types of the blocks and 
some relaxation schemes spanning them. The other schemes 
can be obtained from these figures by rotations. It is conve- 
nient for calculations to subdivide the relaxation processes in 
the block (d) into three parts shown in the dashed box. 

At the stationary state each type of blocks has an ad- 
ditional weight given by the product of the probabilities 
n( h ) of the cells inside the block and the numerical factor 



equal to the numbers of different blocks with the same 
relaxation schemes. Thus, to the blocks on Fig. I the 
following relative weights should be ascribed 

W a = 4n^(n A + n B + n c ) 2 , (9a) 

W b = Anl(n A + n B ), (9b) 

W c = 4nln c , (9c) 

W d = npj. (9d) 

The generating function corresponding to the relaxation 
processes inside the blocks with the weight W a can be 
written with the use of the cr-function, Eq. (||), 

S (N,E,S,W) = 

{a(a(N, 1, 1, W), 1, S, W) - <r(0, 1, S, W) + c.p.}/N a , (10) 

where the notation c.p. implies all possible cyclic permu- 
tations of the symbols N, E, S and W and N a is the normal- 
ization factor such that S Q (1, 1, 1, 1) = 1. To write this 
function we start from the left down cell and define the ar- 
guments of the cr-function corresponding to the toppling 
of this cell. For the process definitely spans the block, 
the left up cell should topple and we write another cr- 
function instead of the symbol N. By going eastward the 
process will terminate inside the block and this branch of 
the toppling process cannot affect the neighboring blocks. 
Hence, we should write the number 1 instead of the sym- 
bol E. The other symbols S and W correspond to the 
branches of the toppling process which immediately get 
out the initial block of cells. 

Analogously, we can write the generating function of 
the blocks with the weight Wb 

S fc (N,E,S,W) = 
{a(cr(N, cr(N, E, 1, 1), 1, W), 1, S, W) - cr(0, 1, S, W) 
+ cr(cr(N,l,l,W),cr(l,E,S,l),S,W) -cr(0,0,S,W) (11) 
+ cr(l, a(cr(N, E, 1, 1), E, S, 1), S, W) - 0, S, W) 
+ C.p.}/N b . 

The generating functions E c and Y,d for the blocks with 
the weights W c and Wd in Fig. [l] are quite similar in 
principle but the expressions are rather long and will be 
published in a complete paper. 

All these generating functions are polynomials describ- 
ing in detail toppling processes inside the blocks. As it is 
easy to check directly, besides the terms with the first 
powers of the symbols N, E, S and W these generating 
functions have also the terms with their squares. The 
first terms correspond to the processes when after the 
relaxation of the block only one particle goes in a given 
direction, while the last ones describe the processes with 
two particles going in the same direction. According to 
the RG ideology, we should consider both of them as the 
transfer of the only new " coarse grained particle" . To 
this end, we can simply replace all the squares of the 
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symbols (n, E, S, w) by their first powers. Finally, we will 
get polynomials of the same form as the original one, Eq. 
(0) , but with the new branching probabilities. 

Now, to write the complete generating function corre- 
sponding to the relaxation of the block 2 x 2 of the cells 
on the lattice £{,, we should average the ^-functions with 
the corresponding weights of blocks defined by the height 
probabilities at the scale b, Eqs. 
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S(N, E, S, W) = {W^ a S a (N, E, S, W) + W 6 S 6 (N, E, S, W) 
+ IF C £ C (N, E, S, W) + W d E d (N, E, S, W)}/N. (12) 

To finish the whole RG scheme we have to define an 
analogous renormalization equations for the height prob- 
abilities n( 2b ) . Due to the nonlocal properties of the dy- 
namics of the sandpile model, the direct renormalization 
of these quantities seems to be very difficult. Instead, we 
will use the stationary state relations, Eqs. j6a[-pd|), to 
define the renormalized height probabilities. 

Given this RG transformation we can study how the 
system evolve under the successive doubling of the scale. 
The final result is that independent of the vectors 
and at the minimal scale the system flows at the 
same nontrivial fixed point. This fixed point is attractive 
in the whole phase space and the system evolves spon- 
taneously toward the critical values p* and n* shown in 
Table | 

TABLE I. Fixed point parameters of the RG transforma- 
tion in comparison with analogous exact results. 
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in the avalanche equal to s varies as P(s) ~ s~ T , for large 
s, where r = 5/4. 

As it has been noted by Pietronero et al. the 
avalanche exponent r can be obtained directly from the 
fixed point parameters. Below, we briefly repeat their 
arguments. By using the discrete length scale = 2 k 
and the avalanche distribution in the form P(r)dr w 
r' 1-2r )<ir we can define the probability that the relax- 
ation process spans the cell of size b^ and dies at the 
neighboring cells not extending over the scale &( fc+1 ) 



K = 



I P{r)dr/ / P(r)dr = 1-2 



2(l-r) 



(13) 



Asymptotically (k — * oo) we can express K in terms of 
fixed point parameters in the following way: 



K=pl(l-n* D )+p* 2 (l-n^) 2 
+ pl(l~n* D f+pl(l-n* D ) i . 



(14) 



Using these two expressions, Eqs. (|l3|,[l4]), the exponent 
t is given by the formula 



T = 1 - 



1 ln(l - K) 

2 m~2 



= 1.248 



(15) 



in excellent agreement with the proposed value r = 5/4. 
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These results obtained with the use of the RG approach 
can be compared with the exact ones for the sandpile 
model. As it has been shown in the paper 0, an 
avalanche can be considered as a sequence of waves of 
topplings each consisting of sites that toppled only once 
in that wave. Being more simple objects, waves admit 
a representation in terms of spanning trees covering the 
lattice sites. The branching probabilities of the span- 
ning trees have been calculated exactly by Manna et al. 
H . Their results are presented in Table |. The hypothe- 
sis that branching probabilities for spanning trees should 
coincide with that for the toppling process seems to be 
quite plausible, but it has yet to be proved. The exact 
height probabilities for the sandpile model presented in 
Table [| were calculated by Priezzhev M . 

Recently, Priezzhev et al. ]To| ] have used the known 
exponents for spanning trees to argue that in the sandpile 
model the probability for the number of sites to involve 
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